library(pacman)
pacman::p_load(tidyverse,data.table,broom,hrbrthemes,plyr,latex2exp,openxlsx,ggpubr,ggpmisc)
rm(list=ls())
################################################################################

###Emissions intensities
df = setDT(read.xlsx('aux_params/parameters.xlsx',  sheet = "emissions_total"))
df = df[, list(spatial_id_name,beef_total_flow_per_t)]
colnames(df)<-c('spatial_id_name','E') 
df_emissions=df

###Change in output 
for(cf_spec in c('dt','certr','certc')){

  if(cf_spec %in% c('dt') ){plot_title='A. No targeting'}
  if(cf_spec %in% c('certr') ){plot_title='B. Region-level targeting'}
  if(cf_spec %in% c('certc') ){plot_title='C. County-level targeting'}
  
  df = setDT(read.csv(file = paste0('results_ic/lf/Q.csv'), header=FALSE))[,list(V1)]
  setnames(df, old=c('V1'), new=c('Q') )
  df$spatial_id_name = df_emissions$spatial_id_name
  df_lf = df
  
  df = setDT(read.csv(file = paste0('results_ic/',cf_spec,'/Q.csv'), header=FALSE))[,list(V1)]
  setnames(df, old=c('V1'), new=c('Q') )
  df$spatial_id_name = df_emissions$spatial_id_name
  df_cf = df
  
  df = merge(df_cf, df_lf, by='spatial_id_name', suffixes = c('_cf',''), sort=F)
  df = merge(df, df_emissions, by='spatial_id_name', sort=F)
  df=setDT(df)
  
  df$x_var=df$E
  df$y_var=100*(df$Q_cf - df$Q)/df$Q
  ylabel=TeX('$\\Delta$ production (pp)')
  
  x_min=min(df$x_var)
  x_max=max(df$x_var)
  y_min=min(df$y_var)
  y_max=max(df$y_var)
  
  x_slope_label = 0.15*(x_max-x_min)+x_min
  y_slope_label = y_min+0.3*(y_max-y_min)
  if(cf_spec=='dt'){x_slope_label = 40; y_slope_label = -67}
  if(cf_spec=='certr'){x_slope_label = 50 ; y_slope_label = -97}
  if(cf_spec=='certc'){x_slope_label = 50; y_slope_label = -143}
  
  spec_m = lm(y_var ~ x_var,  data=df)
  slope_m = round(summary(spec_m)$coefficients[2],2)

  fig_county = ggplot(data = df, aes(x = x_var, y=y_var) ) +
    stat_summary_bin(fun='mean', color='gray40', geom='point', size=0.75)+
    stat_poly_line(formula = y ~ x, se=F, linewidth=0.8, color='black')+
    theme_minimal()+
    labs(title=plot_title, x=TeX('emissions intensity'), y=ylabel)+
    annotate("text", x = x_slope_label, y = y_slope_label, label = paste0("Slope = ", slope_m), hjust=0, size=2.75, color='black')+
    theme(axis.title.x = element_text(size=7, hjust = 1), 
          axis.title.y = element_text(size=7, hjust = 0.95),
          plot.title = element_text(size = 8.25, hjust = 0.5, vjust=2),
          axis.text.x = element_text(size=6),
          axis.text.y = element_text(size=6),
          aspect.ratio = 1,
          panel.grid.minor = element_blank())+
    guides(size = "none")
  
  if(cf_spec=='dt'){fig_mb=fig_county}
  if(cf_spec=='certr'){fig_mb_region=fig_county}
  if(cf_spec=='certc'){fig_mb_targeted=fig_county}

}

figure_final = ggarrange(fig_mb,fig_mb_region,fig_mb_targeted, nrow=1, ncol=3)
ggsave(paste0("../../output/figures/figure5.pdf"), width = 5.5, height = 1.9, units='in')
ggsave(paste0("../../output/figures/figure5.eps"), width = 5.5, height = 1.9, units='in')
ggsave(paste0("../../output/figures/figure5_bw.pdf"), width = 5.5, height = 1.9, units='in')
ggsave(paste0("../../output/figures/figure5_bw.eps"),  width = 5.5, height = 1.9, units='in')

